
N = 101;% c = rand(N,1);     % c_1 + c_2 x + ...+ c_N x^{N-1}
t=linspace(-5,5,N);
c = hann(N).*sinc(t)';

% Roots of Unity
mm = 250;                   % nodes in lower half plane
if mm < N/2+1, error('require mm >= N/2 +1')
else M = 2*(mm-1);
end

rr = [1.1 1 0.9];           % radii of three concentric discs
exps = (0:N-1)';            % exponential factors

fftA = fft(c.*(rr(1).^exps),M); % evaluate polynomial on the rings
fftB = fft(c.*(rr(2).^exps),M);
fftC = fft(c.*(rr(3).^exps),M);

crA = abs([fftA(end) ; fftA(1:mm+1)]); % append end terms and take abs val
crB = abs([fftB(end) ; fftB(1:mm+1)]);
crC = abs([fftC(end) ; fftC(1:mm+1)]);

k = 2:mm+1;                     % indices of middle-disc without end values
Bk = crB(k);                    % value function on the middle-ring

% locate the local minima
locs = (Bk < crA(k-1)) & (Bk < crA(k)) & (Bk < crA(k+1)) ...
            & (Bk < crC(k-1)) & (Bk < crC(k)) & (Bk < crC(k+1)) ...
                & (Bk < crB(k-1)) & (Bk < crB(k+1));
            
indx = find(locs);              % indices of local minima
lx = length(indx);              % number of minima found

lf_rts=[];
if (lx)
    theta = -2*pi/M*(indx-1);   % convert indices into complex numbers
    lf_rts = exp(1i*theta);
end

actual_rts = roots(flipud(c)); % compute roots with companion matrix
plot(lf_rts,'*'), hold on % plot LF approximations,
plot(conj(lf_rts),'*') % companion matrix roots
plot(actual_rts,'or') % and the three concentric discs
uc = exp(1i*linspace(0,2*pi,200)');
plot(uc*rr,'k'), hold off
axis equal, axis off